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Abstract. We report theoretical and simulation studies of phase coexistence in model 
globular protein solutions, based on short-range, central, pair potential representations 
of the interaction among macro-particles. After reviewing our previous investigations 
of hard-core Yukawa and generalised Lennard-Jones potentials, we report more recent 
results obtained within a DLVO-like description of lysozyme solutions in water and 
added salt. We show that a one-parameter fit of this model based on Static Light 
Scattering and Self-Interaction Chromatography data in the dilute protein regime, 
yields demixing and crystallization curves in good agreement with experimental 
protein-rich — protein-poor and solubility envelopes. The dependence of cloud and 
solubility points temperature of the model on the ionic strength is also investigated. 
Our findings highlight the minimal assumptions on the properties of the microscopic 
interaction sufficient for a satisfactory reproduction of the phase diagram topology of 
globular protein solutions. 



1. Introduction 

Phase coexistence in aqueous solutions of globular proteins, like for instance lysozyme 
in water and added salt, has been intensively investigated in most recent years from 
experimental and theoretical points of view (see e.g. [1-7] and references). Such an 
interest is primarily intended to clarify the mechanism of protein crystallization, a 
process whose control, though of crucial importance for the study of protein structure, 
has not yet been characterised in terms of rigorous and predictive experimental 
protocols [8,9]. 

Indeed, the literature documents a diffuse resort to trial-and-error procedures in 
the attempt to grow "good" crystals from the mother solutions [9]; this is also due 
to a lack of an accurate microscopic description of the overall phase behaviour, which 
depends on many factors: for instance, even assuming that a satisfactory description 
of the interactions among "bare" proteins were available, understanding how these 
interactions are altered by the solution variables is a difficult task, also considering 
that specific salt effects determine a ranking of various ions, which is well-known as 
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the Hofmeister series [10]. Despite the strong dependence on solution variables, several 
globular proteins solutions share important features, when crystallization conditions are 
approached. In particular, the second virial coefficient B2 of many protein solutions falls 
in a narrow window of small negative values, whenever the transition takes place [11]. 
In close relationship with such an occurrence, the phase diagram is characterised by 
a metastable protein-protein demixing region located just below the solubility line [4], 
a peculiar property of crucial interest, since experimental and numerical studies have 
shown that the proximity of the critical point to the solubility line should favour the 
formation of good crystals [4, 12, 13]. Besides the crystallization properties, the protein 
demixing binodal plays an important role also in several human diseases [14]. 

These experimental evidences have been rationalized in terms of effective protein- 
protein interactions described through short-range, centrosymmetric potentials [2, 3, 
6,15-17], and such models have been intensively investigated by means of simulation 
and liquid state theories. On the other hand, the rich phase behaviour of short-range 
attractive models has challenged a straightforward generalisation to such systems of 
concepts and methods derived from the theory of simple fluids. In fact, it has been 
recognised (see [18] and references) that the onset of freezing, and generally the overall 
phase behaviour, is substantially affected by the "perturbative" part (with respect to 
the repulsive core) of the potential, rather than being dominated by cxcluded-volume 
and packing (i.e. by entropic) effects, as is the case in the van der Waals picture of 
simple liquids [19]. The observation of two glassy states, originating from different 
cage mechanisms which drive the structural arrest of the system (see [20, 21] and 
references), recently conffi^med by experiments [22], is another aspect of the unusual 
physical properties of systems interacting through short-range forces. 

In this context, one of the most extensively studied model is represented by the 
hard-core Yukawa fluid (HCYF). In fact, early simulations [23] for this system showed 
that, when the decay of the attractive tail is sufficiently fast, the phase diagram is 
topologically similar to that of most crystallizing protein solutions, i.e. with a metastable 
liquid-vapour binodal lying beneath the sublimation line. Yukawa fluids have been 
studied by other authors (sec [24] and references) and by us [25-27], in terms of 
integral equation theories (lET) of liquids [28]. The utility of theoretical approaches, 
along with computer simulation investigations, arises from several circumstances. In 
fact, simulations of systems characterised by short-range attractive forces, as in the 
case of charge-stabilised protein solutions, might be affected by ergodicity problems; 
similarly, the extension to more realistic multi-component cases (where proteins, solute 
ions and water molecules are explicitly considered) faces with severe difficulties, related 
to the strong size asymmetry of the particle species, and to the usually high dilution 
of the macromolccules. Such studies would greatly benefit of the availability of reliable 
theoretical approaches to deal, for instance, with the multicomponcnt configuration of 
the system. We have devoted a considerable effort [25] to test refined lETs against 
rigorous simulation data, starting from the simplest one component HCYF case. 
We succeeded in particular in the optimisation of the thermodynamic consistency, a 
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constraint which is often imposed in the solution of lETs in order to improve their 
accuracy [29]. In this case, we have investigated both the HCYF and a simple DLVO 
model [30] of proteins, widely used in colloid physics. 

As for the determination of the solid-fluid coexistence, we have calculated the 
freezing line of envisaged models [25] in terms of one-phase freezing criteria proposed 
by various authors [31, 32]. We have assessed the accuracy of such criteria in 
reference [33], where theoretical predictions for the Girifalco interaction potential of 
Ceo fuUerene [34], characterised by a short-range decay, have been compared with the 
solid-fluid equilibrium determined from free energy Monte Carlo simulations. Interest to 
such calculations stems from the observation [23] that the Boyle temperature of the Ceo 
model [34] is reproducible with a somewhat short-range HCYF potential. Moreover, the 
Girifalco potential is currently employed as a robust benchmark to test the performances 
of new simulation strategies and theoretical approaches, when applied to fluid with 
short-range interactions (see [35-37]). We summarise in section 2 our most recent 
investigations of the HCYF and the Girifalco potentials. 

We also review in section 2 our results [38] on the crystaUization processes occurring 
in another model for protein solutions, namely the generalised Lennard- Jones potential 
early proposed by ten Wolde and Frenkel [12]. The phase diagram of the model was 
calculated in reference [12] through computer simulation estimates of the free energy, 
and the most favourable conditions for crystal growth from the solution were idcntiflcd. 
We have carried out extended molecular dynamics calculations, in order to compare our 
results for the crystallization kinetics with the free energy calculations reported in [12], 
and with experimental evidences in real protein solutions. 

In section 3 we report on recent investigations [39-41] of the phase diagram of 
lysozyme solutions in water and sodium chloride added salt, modelled in terms of a 
flexible representation of protein-protein interactions, where the forces acting among 
macromolecules are described by means of a DLVO-hke potential [30]. We have been 
prompted to such an investigation since other models — including the HCYF and 
the generalised Lennard- Jones potentials — have hitherto yielded only a qualitative 
reproduction of the experimental phase diagram, and allowed at best to fit the demixing 
line [3,15]; these models also fail to capture the sensitivity of the metastable demixing 
line to solution conditions [6] . We have performed both lET and simulation calculations 
for the protein-rich — protein-poor coexistence curve (the metastable "liquid-vapour" 
binodal), and determined the solubility fine through a simplified version [42-45] of 
the cell theory [46] to estimate the free energy of the sofid phase. We also show 
theoretical results for the dependence of the cloud and solubility temperatures on the 
ionic strength [40], in comparison with experimental data for lysozyme solutions. Our 
conclusions follow in section 4. 
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Figure 1. Model potentials presented in this work. Full line: Yukawa potential 
with z = 7, see equation (Q; dashed line: Girifalco potential for Cgo, equation jSJl; 
dotted line: modified Lennard- Jones potential with a — 50, equation ||2Jl; dot-dashed 
line: DLVO model with a = 3.6 nm, Q = 10 e, Ah — S.GlfceTljo and 6 — 0.18 nm, 
see equations Q-®. In order to compare among different models, all potentials are 
drawn in their reduced units of length and energy, so that v{r/a = 1) = for the 
Girifalco and the modified Lennard-Jones potentials, respectively, and the minimum 
of the well depth is —1. 



2. Theoretical and simulation studies of model fluids interacting through 
short-range forces 

In the last few years we have undertaken an extensive investigation of the phase 
behaviour of the HCYF [25-27], and assessed at the same time the performances of most 
advanced liquid state theories presently available (both semi-analytic and numerically 
solvable). We recall that the hard-core Yukawa potential is written as: 

. X I r < a 

\ —aeexp[—z{r — a)\/r r>a 

where a is the hard-core diameter, e is the potential depth at closest contact and z is 
the inverse screening length (see figure Q). The properties of the HCYF are calculated 
in the context of several liquid state theories, like for instance the modified hypernetted 
chain approximation (MHNC, [47]), the generalised mean spherical approximation 
(GMSA, [49]) and the self consistent Ornstein-Zernike approximation (SCOZA, [48]). 
We shortly recall the basic features of such theories in the Appendix, and refer to 
other papers (see [24, 47-49] and references) for a detailed presentation of their solution 
schemes. 

The theoretical results reviewed here concern an interaction range of the HCYF 
corresponding to a realistic model of colloidal suspensions and protein solutions [2,3,23]. 
We compare our predictions with computer simulations carried out by other authors [23, 
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Figure 2. Phase diagram of the HCYF; temperatures are reported in units of e/Zcb, 
densities in units of a^^, (3 = e/k-oT. Left: z — 7. Liquid-gas coexistence: horizontal 
bars, simulation results [23]; dotted line, SCOZA; full circles, GMSA; open circles, 
MHNC. As=0 locus: diamonds, SCOZA; pluses, GMSA; crosses MHNC; line with 
squares, sublimation line of reference [23]. Right: z — Triangles: GMSA predictions; 
dotted line with diamond, SCOZA binodal and corresponding critical point [51]; line 
with squares: simulation results [23]. See reference [25] for details. 

48,50] or by ourselves. Our investigations [25] show that the MHNC energies are quite 
accurate for z ranging from four to nine; the GMSA energies are reasonably good at 
z = A and show a ~ 10% maximum discrepancy from Monte Carlo (MC) data at the 
highest z. The MHNC and GMSA equations of state and compressibility are also within 
few percents of the simulation results, at not too low temperatures. The SCOZA results 
for the equation of state are quantitatively accurate in any case. At lower z's all theories 
become almost quantitative. 

Results for the phase diagrams are reproduced in figure |2l in the temperature- 
density {T-p) representation, with the freezing lines estimated through the one-phase 
criterion of reference [32], amounting to the fulfilment of the condition that the so- 
called residual multi-particle entropy As vanishes at the phase transition. The limited 
applicability of such structural indicators in the context of "energetic" fluids [18] has 
been carefully discussed in reference [33]. As visible in figure at ^ = 7 the GMSA 
qualitatively reproduces the relative location of the binodal and sublimation lines 
estimated in simulations of reference [23]. SCOZA results look fairly good; such an 
accuracy seems useful for the investigation of the z = 9 binodal, where no simulation 
results are available. Also in this case the SCOZA binodal [51] lies beneath the GMSA 
sublimation line; the latter one in turn, fairly reproduces the high density portion of the 
sublimation line. It appears that the combination of GMSA and SCOZA may provide 
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Figure 3. Helmholtz free energy as a function of the system volume along supercritical 
isotherms for the HCYF with z = 7 at ksT/e = 0.5 (left), and a DLVO model with 
tj = 3.6nm, Q = lOe, Ah = 8/cb720 and 6 = 0.144nm [see equations for 
symbols] at T = 442.5 K (right). Open circles: MC simulations; dotted hne: MHNC; 
fuU hne: HMSA. See reference [29] for details. 



a description of the relative position of phase coexistence hnes of the HCYF in the very 
short-range interaction regime. 

In a subsequent study [29], we have exploited the sensitivity of integral equation 
approaches like the MHNC and the Hybrid Mean Spherical Approximation (HMSA, [52], 
see the Appendix), to the procedure adopted in order to achieve the (partial) 
thermodynamic consistency of the theory. We have analysed in particular the 
application of MHNC and HMSA to the one-component HCYF and to a system 
interacting through a short-range DLVO potential [30] (see figure ^ and the next section 
for details). We have found that the enforcement of a global consistency constraint, i.e. 
a procedure in which one imposes the equality of the equation of state, as calculated 
along different routes from structure to thermodynamics, is in general more accurate 
than the requirement of a local consistency, based on the equality of the isothermal 
compressibilities. The improvement associated to the global consistency becomes crucial 
in the short-range interaction regime; in particular, as visible in figure El the Helmholtz 
free energy and the chemical potential are almost quantitatively predicted in the global 
consistency approach [29]. These results may help in the selection of integral equation 
approaches appropriate to determine in a confident manner the phase diagram of both 
pure fluids and binary mixtures of model short-range systems. 

As far as a theoretical characterisation of solid-fluid equilibrium is concerned, we 
have recently investigated [53] the solid phase of a pair, short-range model introduced 
by Girifalco [34] to describe particle interactions in Ceo: 
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Figure 4. Fluid-solid coexistence densities (left) and pressures (right) of the Girifalco 
Ceo model. PT predictions are used for the sohd phase and MC data for the fluid 
phase [33]. Open symbols refer to different prescriptions for the effective hard-core 
diameter entering the theory — see equation (refeq:sigmal) in the Appendix. Squares 
are full Monte Carlo results [33]. The Gibbs Ensemble Monte Carlo liquid- vapour 
coexistence line [57] is also shown. Dashed lines are guides to the eye. See reference [53] 
for details. 



where s = r/d, ai = N'^A/12d^, and a2 = N'^B /90d^'^; N and d are the number 
of carbon atoms and the diameter, respectively, of the fuUerene particles, A = 32 x 
10~^° ergcm^ and B = 55.77 x 10~^°^ ergcm^^ are constants entering the Lennard- 
Jones 12-6 potential through which two carbon sites on different spherical molecules 
are assumed to interact. For Ceo, d = 0.71 nm while the node of the potential (j2)), 
the minimum, and its position, are tq ~ 0.959 nm, e ~ 0.444 x 10^^^ erg and 
'^min = 1-005 nm, respectively (see figure [T]). The phase diagram of the model has 
been recently characterised in terms of MC calculations of the free energies of the fluid 
and solid phases [33,54]. We have also determined theoretically the free energy of the 
solid phase in terms of a standard perturbation theory (PT, [28], see the Appendix). 
As visible in figure |31 such an approach yields reliable predictions for the coexistence 
properties and, more generally, for the whole solid phase behaviour. Since the PT 
embodies a first-order expansion of the model's free energy around the free energy 
of a reference hard-sphere crystal, our results noticeably imply that the hard-sphere 
properties dominate, to a large extent, the structure of the solid phase; conversely, 
as discussed in [33], the structure of the fluid phase is markedly affected by energetic 
aspects, related to the attractive part of the potential. Our results further support the 
use of PT to characterise the solid phase of systems interacting through short-range 
forces, already documented in references [21,55,56]. 

We have also investigated in a recent study [38] the phase transformations occurring 
in a system of spherical particles interacting through a short-range potential, obtained 
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Figure 5. Generalised Lennard- Jones potential — see equation (PJ- Snapshots of MD 
configurations at progressively decreasing temperatures (in units of t/k-B), together 
with the corresponding radial distribution functions; (a) nietastable homogeneous fluid; 
(b) low-density bubble in a high-density fluid, just below the metastable liquid-vapour 
separation; (c) crystal nucleus in the high-density fluid; (d) extended defective crystal 
phase. See reference [38] for details. 
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as a generalisation of the Lennard- Jones interaction law [12]: 
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By changing the parameter a, this expression gives rise to a family of potentials whose 
shape ranges from a Lennard- Jones like potential (for a ~ 1) to potentials having a steep 
repulsive part, followed by a deep and narrow attractive well, approaching the sticky 
sphere limit with increasing values of a. In our computations, we set a = 50, thus 
recovering the model early proposed in reference [12] to study the phase behaviour of 
globular protein solutions (see figureH]). We have focused on the crystallization processes 
of this system [38]. In particular, we have performed extensive molecular dynamics 
(MD) simulations of the metastability induced by progressively cooling the fluid below 
the fluid-solid and the (metastable) liquid-vapour coexistence at fixed density. We 
have analysed in detail how the crystallization process evolves in the simulated system, 
and the interplay of crystallization with the metastable liquid-vapour separation. In 
figure El we show a sequence of snapshots of MD simulations at pa^ = 0.5 (i.e. a density 
moderately higher than the critical density of the liquid- vapour binodal), for various 
decreasing temperature in a path which crosses the binodal line at a reduced temperature 
T ~ 0.41. Model (jH)) for globular proteins is admittedly a crude one. Nonetheless, 
qualitative similarities with the experimental behaviour still survive, as for instance the 
enhancement of the crystallization kinetics due to the presence of the metastable liquid- 
vapour separation; this phenomenon is very likely determined by the presence of the 
high density fluid, which decreases both the interfacial free energy required to nucleate 
the crystal, and the amplitude of the density fluctuations, necessary to reach the solid 
density. Another interesting outcome is that the simulated crystallization process does 
not appear to be accompanied by any precursor effect: in fact, when any of the criteria 
we devised to this aim was met, the nucleation process was in any case already well 
underway. Our observations do not strictly prove, however, that precursors do not 
exist, but only that the chain of events leading to the transition involves the non-trivial 
coupling of several variables, as shown experimentally in reference [58]. 



3. Colloidal models for lysozyme solutions 

In a recent series of studies [39-41] we have considered a prototype globular protein 
solution, namely lysozyme in water and sodium chloride added salt, extensively 
characterised from the experimental point of view (see e.g. [1-7] and references). In 
particular, we have parameterised the physical properties of the protein solution in the 
framework provided by the DLVO theory of charged colloidal suspensions [30]. We recall 
that the DLVO potential of mean force is given by the sum of a repulsive, Debye-Hiickel 
contribution. 



Phase coexistence in model protein solutions 



10 




Figure 6. Interaction factor as a function of the solution ionic strength for lysozyme 
in water and NaCl added salt. Circles: SIC [7]; squares: SLS [2,6]; diamonds: SLS [1]; 
triangles up: SLS [5] at T = 303.15 K; triangles down: SLS [5] at T = 293.15 K. Full 
line: DLVO potential with An = 8.61fcB'72o. See reference [41] for details. 



and an attractive part, represented by a short-range van der Waals term: 

An 



vuA[r) = 



(j^ (T^ r'^ — (J^ 



(5) 



The resulting total interaction is: 



f oo r < a + 5 

VDLYo{r) , ■ (6) 

[vuA{r) + VBuir) r>a + d 

Here a represents the effective diameter, Q = ZpC is the net charge (in electron units) 
of the macroparticle and A-^ is the Hamaker constant; and eo are, respectively, 
the (solution) relative and the vacuum dielectric constants; xdh is the inverse Debye 
screening length. The Stern layer thickness S, related to the intrinsic size of counterions 
which condense on the protein surface, is introduced in equation © to circumvent the 
singularity of the attractive term at r = a [1]. 

In our analysis, we have fixed all parameters entering equation ©, but the Hamaker 
constant, to values commonly accepted on the basis of several experimental evidences 
for lysozyme i.e. a = 3.6 nm, Q = lOe, and 6 = 0.18 nm [1,59,60]. The Hamaker 
constant is used as a free parameter to fit the experimental interaction factor defined 
as: 

POO 

Ks = Ki"" + 24 / x^l- exp[-/3t;DLVo(a:)]} dx , (7) 

where x = {r/a — 1) and Kg^^ = 8 is the interaction factor corresponding to a fluid 
of hard spheres. The interaction factor is proportional to the second virial coefficient, 
B2 = Ksi'/2M, where u = 0.703 ml/g and M = 14300 Da are respectively the specific 
volume and the atomic mass of lysozyme. As shown in figure El we best-fit Self 
Interaction Chromatography (SIC) results, reported in reference [7], as well as Static 
Light Scattering data from various sources (SLS [1,2,5,6]), for a lysozyme solution in 
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the range [0.1 — 0.6] M NaCl at a pH close to that of the experimental phase diagram [4]. 
Our fits are carried out in a temperature range relatively far from the demixing region; 
under these conditions the second virial coefficient of the lysozyme solution is expected to 
vary moderately, and almost linearly, with the increase of the temperature. An average 
Hamaker constant is then obtained. An = 8.61 (in ^6^20 units) with a maximum half- 
dispersion equal to 0.26. 

The extension of a "colloid physics approach", as the DLVO theory, to globular 
protein solutions has long been debated. In fact, the literature reports either successful 
applications (see e.g. [17,60-62]), or evidences of opposite sign, concerning the accuracy 
and the theoretical foundation of this picture [63-65]. Our observation is that, in 
agreement with previous studies, the validation of the present approach mainly derives 
from a positive, non-trivial, parameterisation of phenomenological observations. In 
particular, we envisage the functional form of the attractive potential in equation ^ 
as flexible enough to accommodate not only the van der Waals forces, but also other 
non-DLVO interactions, as hydration effects [60,61,66]. We assume in fact that these 
forces act to artificially infiate the value of the Hamaker constant, which plays a different 
role from the original one, only related to dipolar attractions. Within this framework, 
specific effects, related for instance to the salt identity, are taken into account only 
indirectly, through the fit of experimental data. Although one could introduce additional 
interaction terms (which can be salt-specific), in order to increase the fiexibility of the 
phenomenological potential, we shall show that a fine tuning of the van der Waals 
amplitude is sufficient to set the proper order of magnitude of the effective protein- 
protein interaction. In particular, this attractive term is the leading contribution to the 
effective potential because repulsive electrostatic interactions are almost screened out at 
biological salt molarities, thus constituting a minor correction to the potential of mean 
force. 

On the basis of the above model, we have calculated the phase diagrams of the 
lysozyme solution at Ig = 0.51 M reported in figure [7| The protein-rich — protein-poor 
phase coexistence line is determined by means of Gibbs Ensemble Monte Carlo [67] 
simulations, carried out on a system composed of 1024 particles (with a 3a interaction 
cutoff). As for the solubility line, the free energy of the fluid phase is calculated in 
the framework provided by the reflned HMSA theory [52] whereas, for the solid phase, 
we have used for the chemical potential a simplifled expression [42-45] derived in the 
framework of the cell theory [46], 



This equation provides a direct link with the essential properties of the protein crystals, 
namely the average number of contacts, Ug, and the translational freedom along one 
axis, A [with A* = \/{a + 6)], of the protein inside the unit cell, so that A'^ constitutes a 
rough estimate of the volume accessible to the protein inside the cage formed by its flrst 
neighbours; fiQ is the standard part of the chemical potential, and eEFF is the minimum of 
the effective potential ©. A constant value log[l^/(cr + 5)^] [45] has been included in /xq. 




(8) 
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Figure 7. Phase diagram of a lysozyme solution in water and NaCl at Ig = 0.51 M. 
Open circles and full triangles: cloud and solubility points, respectively [4]. Error 
bars with diamond: DLVO demixing envelope and corresponding critical point. Full 
and dashed lines: DLVO crystallization curves with parameters for equation © in the 
legend. 

where Vp is the protein volume. The coexisting fluid branch is then found by equating 
the chemical potentials of both phases [43]. We fix the number of contacts between 
nine and ten, in agreement with the estimates of reference [68] ; a comparable number of 
contacts has been specifically reported for lysozyme crystals in [69]. Our results for the 
phase diagram are compared in figure [7| with the experimental demixing and solubility 
lines reported in reference [4]. As visible, the theoretical demixing envelope coincides 
with the experimental locus of cloud points with a critical temperature T^r — 284 K, 
very close to the experimental outcome Tcr = 282.5 K. The theoretical crystallization 
boundary, calculated with = 9 and A = 0.03 nm in equation (jH)), is also in good 
agreement with the experimental curve. We obtain similar results if another plausible 
number of contacts, namely Us = 10 (with A = 0.017 nm) is assumed (see figure Hj). 
All values of A are compatible with the temperature factors of the Bragg intensities in 
lysozyme crystals, where the typical range for A is [0.014 — 0.033] nm [70]. 

As visible from figure |H1 satisfactory results are also obtained for the variation of 
cloud and solubihty temperatures as a function of the ionic strength. These results 
are reported in a very recent work [40] where we have fitted the Hamaker constant 
against the collective diffusion coefficient measurements carried out by Beretta and 
coworkers [60]. It appears that the theory substantially capture the experimental 
dependence [4, 10] of the critical and solubility temperatures on the logarithm of the 
ionic strength. 

In summary, although a protein model where we only consider central forces might 
appear oversimplified, a reasonable prediction of the phase diagram of the lysozyme 
solution may be obtained through a "colloidal" model, provided a judicious fit of the 
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Figure 8. Critical (left) and solubility (right) points as a function of the ionic 
strength for lysozyme solutions at fixed protein concentration p = 90 g/1. Full symbols; 
experimental data from reference [4] (squares), [10] (circles, pH=7.8), and [4] (diamond, 
pH=4.5). Open symbols: DLVO results. Curves are linear fits of theoretical and 
experimental data. 



experimental second virial coefficient of the solution [1,2,7] is performed. 

On the other hand, it is evident that a detailed analysis of the complex behaviour 
of protein solutions would require a refined, highly specific treatment of the microscopic 
interactions; several improvements to the DLVO-like model have been recently proposed 
as, for instance, different shape [71] or multi-site representations of protein interactions 
(sec e.g. references [43,72-74]). In comparison with such studies, our results support the 
observation that in the fluid phase, close to the critical point, a distributed approach, 
as envisaged for instance in [72] may be well approximated by its spherically averaged 
counterpart. 

4. Conclusions 

We have reported results of extensive investigations of short-range potential models 
suited to describe, in a realistic albeit qualitative manner, the physical properties of 
prototype globular protein solutions like for instance lysozyme in water and added salt. 
We have focused on the phase diagram of these systems, and first reviewed predictions 
from simulations and integral equation theories of the fiuid state for the hard-core 
Yukawa fiuid. We have carried out a systematic assessment of various refined theories — 
like the MHNC, the SCOZA and the GMSA — against simulation data, and investigated 
the improvement of the strategy adopted to enforce their thermodynamic consistency. It 
emerges that integral equations can reach a good degree of accuracy in the prediction of 
phase coexistence properties. In particular, the basic property of the phase diagram in 
many protein solutions, namely a liquid-vapour coexistence line metastable with respect 
to the fiuid-solid equilibrium, is correctly reproduced. Hopefully, these theoretical 
approaches may be used as a complementary tool, especially when simulations are 
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applied with difficulty, as in the treatment of multicomponent systems characterised by 
the high dilution of one of the components, or by strong charge and/or size asymmetries. 

We have also reviewed a recent molecular dynamics investigations of crystallization 
processes in a generalised Lennard- Jones model of protein solutions. This study shows 
that some features of the true crystal formation kinetics are qalitatively represented, 
even if a reahstic description of the natural process (especially as far as the crystallization 
time scale is concerned) is still out of the possibihties offered by such an approach. 

We have then reported on an investigation of protein solutions based on the DLVO 
potential, widely adopted in the representation of colloidal systems. We have shown 
that the phase coexistence lines of lysozyme solutions in water and sodium chloride 
added salt can be reproduced through a combined use of theories and simulations, when 
the DLVO potential is parameterised on the experimental second virial coefficient of the 
protein solution at moderate ionic strengths. Preliminary results for 7-crystallin in water 
and added sodium phosphate also turn out quantitatively accurate. Our procedure is 
therefore potentially suited for a variety of protein systems, and presently its extensive 
test appears to us as essentially hampered by the scarcity of phase diagram and B2 
measurements in the same solution conditions. 

In conclusion, basic assumptions on the effective interaction among globular 
proteins in solution are sufficient to accurately reproduce the phase boundaries of such 
systems. In particular, a parameterisation of a central, pair potential, based on the 
experimental second virial coefficient allows for a confident prediction of the true cloud- 
points and solubility lines. Since a much less effort is required for the experimental 
determination of B2, the advantage of the procedure here exploited becomes evident. 
Our predictions are accurate also for high protein concentrations, although experimental 
measurements of B2 are carried out in the dilute regime. These findings agree with the 
hnearity observed in several Debye plots of scattering data [1], indicating that B2 values 
measured in under-saturated solutions slightly change in the super-saturated regime. 

Appendix: Theoretical methods 

The properties of model fiuids discussed in this work have been determined theoretically 
in the context of integral equation theories of hquids [28] hke the modified hypernetted 
chain approximation (MHNC, [47]), the generalised mean spherical approximation 
(GMSA, [49]), the self consistent Ornstein-Zernike approximation (SCOZA, [48]), and 
the Hybrid Mean Spherical Approximation (HMSA, [52]). In the MHNC scheme one 
considers the Ornstein-Zernike equation 



(where p is the density of the fluid, g{r) is the radial distribution function, rdf, 
h{r) — g{r) — 1 and c(r) are the pair and direct correlation function, respectively) 
coupled to the formally exact expansion result for the rdf 




(A.l) 




(A.2) 
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where v{r) and B{r) are the pair potential and the bridge function, respectively, of 
the system under study, while (3 = l/k^T with the Boltzmann constant and T 
the temperature. A closure to equation (jA.ljl is then obtained through equation ()A.2jl 
by approximating the bridge function for the potential v{r) in terms of that of a hard 
sphere fluid characterized by a hard core diameter a* chosen so as to enforce the (partial) 
thermodynamic consistency of the theory. 

In the GMSA, a solution to the Ornstein-Zernike equation (2) is obtained by hinging 
on the exact core condition g{r) = for r < a, descending from the fact that v{r) = oo 
on such a distance range, coupled to the approximation for the direct correlation function 
outside the hard-core diameter: 

c{r) = —j3v{r) + K exp[—z{r — a)]/r r>a (A. 3) 

where K and z are used as adjustable parameters in order to impose the internal 
thermodynamic consistency of the theory. 
In the SCOZA scheme c(r) is written as: 

c(r) = -A/3v{r) + Khs exp[-2Hs('^ - ^)]/^ r > a , (A.4) 

Kus and zhs are predetermined by setting t>(r) = in equation ()A.4|1 and requiring that 
both the compressibility and the virial routes to thermodynamics yield the Carnahan- 
Starling equation of state for a hard-sphere fluid. The factor A is then obtained by 
imposing that the compressibility and energy routes yield the same thermodynamics. 

The HMSA is so called because it interpolates between the hypernetted chain 
(HNC) [24] approximation and the soft mean spherical approximation (SMSA) [24]: 
it is written as 

In equation ()A.5|) the potential v{r) is splitted into a soft core repulsive part, [fi(r)], 
and an attractive part, [f2(r)], according to: 

viir) = { ^ V2ir) = { ^ ^ (A.6) 

(0 r > rmin [ v{r) r > rmm 

where r^^^ denotes the position of the minimum of potential v{r) and /(r) = 1 — exp[,^r]. 
The parameter ^ is used to enforce the thermodynamic consistency. 

In order to determine theoretically the properties of the solid phase of the Girifalco 
Geo model we have employed a standard perturbation theory (PT, [28]). In this scheme, 
the model's free energy per particle, /, is calculated as a first-order expansion around 
the free energy of a reference hard-sphere solid, /hs, according to: 

/?/ = /5/hs + yJ ^P-t(^) ^Hs(^) dr , (A.7) 

where ^hs(^) radial distribution function of the reference hard-sphere solid. In 

order to determine the perturbative part of the potential, fpert('") and the diameter a* 
of the reference hard-sphere system entering equation ()A.7|) . we have followed the WGA 
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prescription [19], that amounts to split the original potential according to equation ()A.6|1 
with Vpert{r) = V2{r) and fref = Vi{r); then [75]: 



a 



/ {l-exp[-/3t;,ef(r)]}dr. (A.8) 

JO 
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